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Abstract 

In this paper, we investigate a non-iterative imaging algorithm based on the topological 
derivative in order to retrieve the shape of penetrable electromagnetic inclusions when their 
dielectric permittivity and/or magnetic permeability differ from those in the embedding 
(homogeneous) space. The main objective is the imaging of crack-like thin inclusions, but 
the algorithm can be applied to arbitrarily shaped inclusions. For this purpose, we apply 
. multiple time-harmonic frequencies and normalize the topological derivative imaging func- 
o I tion by its maximum value. In order to verify its validity, we apply it for the imaging of 
two-dimensional crack-like thin electromagnetic inhomogeneities completely hidden in a ho- 
mogeneous material. Corresponding numerical simulations with noisy data are performed 
for showing the efficacy of the proposed algorithm. 
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1. Introduction and preliminaries 

The main objective of this paper is the development of a topological derivative based one- 
step iterative imaging algorithm for thin electromagnetic inclusions completely embedded 
in a homogeneous domain, via boundary measurement. For proper beginning, we review 
related mathematical models, and corresponding formulas, followed by a brief condensation 
of recent results and an outline of the current paper. 

Let f2 C be a homogeneous domain with smooth boundary dQ, which is a curve; 
this domain contains a thin, curve-like homogeneous electromagnetic inclusion. Let us as- 
sume that this thin inclusion (denoted as F) is represented in the neighborhood of a simple 
smooth curve a := cr(x) as 

F = {x + 7n(x) : X G a, 7 G (-/i, h)} , 

where n(x) is the unit normal to a at x and /i is a positive constant that denotes the thickness 
of F refer to Figure [H Throughout this paper, we assume that the applied frequency is of 
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Figure 1: Two-dimensional thin electromagnetic inclusion T of thickness 2h. 



the form u = ^ for the given wavelength A, the thickness /i of F is sufficiently small with 
respect to A (/i ^ A), and the inclusion does not touch the boundary dQ so that it must be 
located at some distance from dfl. In other words, there is a nonzero positive constant s 
such that 



Let every material be classified by its dielectric permittivity and magnetic permeability 
at a given frequency u. Let < £o < +00 and < iiq < +00 denote the permittivity and 
permeability of the domain Q, and < e < +00 and < < +00, those of the inclusion 
r. Then, we can define the piecewise constant dielectric permittivity £(x) and magnetic 
permeability /i(x) as 



dist((7, dVl) = s 3> /i. 




for X G r 



for X G n\T 



and 




for X G n\T 



for X G r 



(1) 



respectively. For the sake of simplicity, we set £q = fiQ = 1, e > Eq, and fi > fiQ. 

At a given frequency u, let u^''\'x.;u) be the time-harmonic total field satisfying the 
Helmholtz equation in the existence of F, 




with transmission conditions 



u 




and 



1 <9M(')(x;a;) 



+ 



1 aM«(x;a;) 
/i drjix.) 



on dr. 



Here, and r7(x) represent the unit outward normal to x G dQ and x G dT, respectively, 

2 



subscript ± denotes the limiting values as 
"^(xJli = lim 'u(x±t77(x)) and 



= lim \ , " for X G dT, 



and = (cos sin 5;) denotes a two-dimensional vector on the unit circle Similarly, 
let 'Ubjc(x;a;) = e^'^'^''^ denote a field satisfying ([2]) without F, i.e., a background solution. 
Throughout this paper, we assume that is not an eigenvalue of ([2]). 

As mentioned earlier in this section, the main purpose of this paper is to develop a fast, 
non-iterative algorithm for imaging a thin inclusion F completely embedded in a domain Q, 
via the boundary measurements ^('^(x; w), x G dfl. Note that there is a remarkable number 
of interesting inverse scattering problems for reconstructing thin electromagnetic inclusions 
and/or perfectly conducting cracks hidden in a structure (such as bridges, concrete walls, and 
machine constructions) from boundary measurements, refer to js], 0] and references therein. 
For this purpose, various iterative and non-iterative imaging algorithms have been developed 
and successfully applied to various problems, for exam ple, level-set method ji], [Til, 28 



Multiple Signal Classification (MUSIC)-type [61 Isl , lid. l25l l27l| . linear sampling method 



141 Il9| and multi- frequency based algorithms [6|, 118|, |21|, |22|, |26[. From many researches, it 



turns out that non-iterative imaging algorithms are fast, simple, effective, and extendable to 
multiple targets; however, they require a large number of incident directions and boundary 
measurements. In contrast to the non-iterative algorithms, iterative imaging algorithms do 
not require a large number of incident directions and boundary measurements. However, 
they require complex calculation of the so called Frechet derivative, adequate regularization 
terms for each iteration step, a good initial guess whose shape is close to the unknown 
target (here, F) and a priori information of target, e.g., material properties, thickness, 
location. Owing to these considerations, the realization of a trade-off between non-iterative 
and iterative imaging algorithms is an interesting research topic. 

Topological derivative strategy has been developed for this purpose. Recently, this strat- 
egy was successfully applied to the shape optimization and imaging of small and crack-like 



inhomogeneities, see |5|, lid, ll2|, ll3|, lly, |23|, |2J, |3l[ for instance. For successful application of 



this strategy theoretically, a large number of incident directions and corresponding scattered 
fields are required. Unfortunately, for practical application, it is extremely difficult to in- 
crease the number of such fields owing to the high configuration costs, unavoidable random 
noise, and so on. 

The above limitation has motivated us to consider an improved topological derivative for 
imaging thin, extended electromagnetic inclusions. For this purpose, we propose an imaging 
functional based on the topological derivative at multiple frequencies. We explore some 
properties and limitations of traditional topological derivative based imaging functional, 
and we aim to improve them accordingly. 

The remainder of this paper is organized as follows. In section |2l we briefiy introduce 
the topological derivative based imaging functional derived in [2^. A normalized multi- 
frequency imaging functional is proposed in section [31 In section HI we present the results 
of numerical simulations to illustrate the advantages and disadvantages of the proposed 
imaging algorithm. Finally, we conclude this paper in section [51 

3 



2. Review of normalized topological derivative at single frequency 

In this section, we shall introduce the basic concept of topological derivative operated 



oi top e 

at a fixed single frequency. We would like to mention [s], Ho, ll2|, |13|, |16|, |23|, |24, |31| for 
detailed discussions. Let UtJt{^;u) and u^}^{'x.;u) be the total and background solutions of 
dlj), respectively. The problem considered herein is the minimization of the following energy 
functional depending on the solution M*^'^(x;a;): 

1=1 1=1 "^^^ 

(3) 

Assume that an electromagnetic inclusion S of small diameter r is created at a certain 
position z G Q\dQ, and let Q\Ti denote this domain. Since the topology of the entire domain 
has changed, we can consider the corresponding topological derivative dx^lz) based on E,{fl) 
with respect to point z as 

dTE(z;u)= lim ^ ' ' / , , 4 

where (/?(r; u) — )■ as r — )■ 0+. From we can obtain an asymptotic expansion: 

E(^]|S; uj) = E{Q; u) + if{r; u)dTE{z; u) + o{^{r; u)). (5) 

In [iS], the following normalized topological derivative imaging function ETD(z;aj) has 
been introduced: 

if dTEs{z]u) ^ dTEfj_{z]Uj) 



2 \max[o(rllle(z; wjj max[aTE^(,z; wjj / 

Here, (iTE£(z;a;) and dxE^^z^u) satisfying ([5]) for purely dielectric permittivity contrast 
{e 7^ £o and /i = fio) and magnetic permeability contrast {e = Eq and ;U 7^ /io) cases, 
respectively, are explicitly expressed as (see (24J] ) 

dTE,(z; cu) = ^He ^ [^;S(z; u;)Mlil(z; w)") , (7) 
1=1 ^ ^ 

drE^iz; u) = mtJ2 (v^fS(z; ^) " Vwlil(z; a;)") , (8) 



where t;fdj(x;a;) satisfies the adjoint problem 

At;S(x;a;)+u;2^S(x;a;) = in Q 

^|^ = .£)(x;c.)-.a(x;c.) on 9a 

Some remarkable properties of O and ([H]) for small and extended thin electromagnetic 

inclusions can be found in ^ and j^, respectively. 
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3. Introduction to normalized multi-frequency topological derivative: theory 
and calculation 



Topological derivative based imaging algorithm is well known for its fast imaging per- 
formance and robustness with respect to random noise (see [5| for instance). However, 
when measured data is affected by a considerable amount of noise and/or the number of 
incident directions L is small (see 2J] for the effect of L), one cannot obtain a good re- 
sult. In order to address these issues, we refer to multi- frequency based imaging techniques 
[s, 17, 18, 21, 22, 26|, and we consider the following normalized multi-frequency based topo- 
logical derivative imaging function: for several frequencies {uk : k = 1,2, ■ ■ ■ , K}, define 



1 ^ 1 ^ 

E(z; K):=-J2 Eto(z; u^k) = 



k=l 



dTEJz;Uk) 



+ 



dTEJz;Uk 



2K ^ Vmax[(iTEe(z; Uk)] max[dTE^(z; Uk)] J ' 

), respectively, for oj 



(10) 



u)k, k 



where dx^si'^'T^k) and (iTE^(z; Wfc) satisfy ([7]) and 
1,2,---,K. 

From now on, we will analyze the properties of fllUp . For this purpose, we recall the 
following result from Note that only a concise proof of Lemma [XT] is introduced in 24 1: 
we have provided a detailed proof of Lemma 13.11 in Appendix [Bl 

Lemma 3.1. Let A B imply that there exists a constant C such that A = BC , and let 
yic{f) denote the real part of f . Then, ^ and ^ satisfy 



dTEJz;(^k) 



dTEf,{z;Uk 



1=1 



1=1 



e - £oje 
1 1 



/i /io 



di ■ t(x) + 2 



1 



/i 



/io /io 



d; ■ n(x) 



where t(x) and n(x) are unit vectors that are respectively tangent and normal to the sup- 
porting curve aat^. 

With this, we can obtain the following result. 

Theorem 3.2. Assume that the number of incident directions L{> 4) is small and the 
applied number of frequencies F is finite (F < +oo^; then, /[W\) becomes 



E(z; K)^- 



Ei(z;K) 



+ 



E2(z;K) 



2 Vmax|Ei(z;K)| max |E2(z; iiT)]/ ' 
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where 



1=1 •^'^ ^ 

1—1 J cr 



d/ ■ (x — z) ) cos 



■ (x - z) )(ia-(x) 









-^)d,.n(x)^ 


K-- 




)dz-t(x) + 2f— 











X Jo 



/ UK - ^1 

V 2 



d; ■ fx — z) COS 



d; ■ (x — z) J dcr(x 



and jo{x) denotes the spherical Bess el function of order zero, 

sin a; 



Jo (a;) 



Proof. First, we consider the term Ei(z; K) by evaluating 
y^^dTE^{z]Uk) - / dTE^{z)duj ^ i£-£o) 



fe=i 



^ia;di-(x-z) 



(icr(x) 



du 



1=1 



J2 /(e-£o)(^e / e''^'^'-^''-'Uuj]d(T{^). 



(11) 



Performing an elementary calculus yields 



gia;di-(x-z) 

zd; ■ (x — Z 
1 



git^xdr(x-z) _ gia;idr(x-z) 



zd; ■ (x — Z 

1 ' ^ 

cosfw/i'di ■ (x — z)) — cosfwid; ■ (x — z)) 



id; ■ (x — z 

+ i sinfw/^d; ■ (x — z)) — i sinfwidi ■ (x — z)) 



d; ■ (X - Z ) 



COS 



d; ■ (x — z) sin 



{uk - ^l) 



d; ■ fx - z] 



+ i sin 



{uk + ^l) 



d; ■ fx — z) sin 



di • (x - z] 



Therefore, by taking the real part of the above formula, flTT]) can be approximated as 



K ^ r 

k=i 1=1 '^'^ 



sinf^id; ■ (x - z)) 
d; ■ fx - z) 



cos(^2di ■ (x - z))dcr(x) 



2 ^ f 

- ^ 1 {e- eo)jo(6di ■ (x - z)) cosf^d; ■ (x - z))c/a(x), 



(12) 



where 



4i := 7, and 42 : = 



2 2 
Hence, by taking the maximum value of ( [T2|) and using it for normahzation, we can obtain 
the desired structure of ]Ei(z; K). 

Next, we consider the term E2(z; K). Since d^, t(x), and n(x) do not depend on ojk-, we 
can similarly obtain the following approximation 



K 



k=l 



E 



1=1 



2( i - — Id/ -trx) + 2' ^ ^ 



I ■ Li^x; ^\di-n{ 



sin(^id, ■ (x - z)) 



d/ ■ (x - z) 



cos(^2dj ■ (x - z))(ia-(x) 
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E 



Z=l 



21 - --)d,-t(x) + 2( - -4 )drn(x) 



io(6d« ■ (x - z)) cos(^2di ■ (x - z))dcr(x) 

(13) 

By applying the maximum value of f[T5]) . the structure of IE2(z; K) can be obtained. 

Theorem 3.3. Assume that the number of incident directions L is sufficiently large and the 
applied number of frequencies F is finite (F < +00); then, /[TU\) becomes 



2 Vmax|E3(z;K)| max |E4(z; 



with 



E3(z; K) = 27r j^{e - Eq) (^A(t; uk) - A(t; Ui)^ da{^) 











K- 


--')d,-t(x)- 













Here, A(t; u) is defined as 



A{t;u) := uJoicut) + ( JiM)i^oM) - Joi(^t)Hi{ujt) 



(14) 



where Jn{x) denotes the Bessel Junction of order n of the first kind and Hn denotes the 
Struve function of order n (see Jj, Chapter 11]). 



Proof. By employing the result in 17|, Lemma 4.1], the following relation holds: for suffi- 
ciently large L, 



V- gia;dr(x-z) ^ j ei^<^-{^-^)dS{d) = 27r Jo(^|x - z|). 

1=1 



(15) 



Let K — oo; then, applying an indefinite integral of the Bessel function (see [30|, page 3]), 
J Jo{t)dt = tJoit) + y (^Mt)Ho{t) - Jo{t)H^it)^ , 



yields 



K 



K 



{e - eo)Jo{uk\^- z\)da{x.) ^ 2t[ / {e - eo)Jo{u\:si - z\)duda{x) 

= 2t: j{e- Sq) (^k{t- ujk) - A(t; Wi)^ rfor(x), 

where function A(t; a;) is given by (IT^ . Hence, we can obtain the structure of IE3(z; K) via 
the above identity. Similarly, the structure of ]E4(z; K) can be identified. 

Theorem 3.4. Assume that the number of incident directions L and frequency ojk ore suf- 
ficiently large enough, and K is infinite (K — )■ oo); then, [TP) becomes 



E(z; K) 



where 



E5(z; K) = !{e- e^)-^^^da{^) 

Ja |X-Z| 



2 \max |E5(z; max |E6(z; A' 



27r 



E6(z;K) 



1 1 



2 di-t(x) + 2 ^ di-n(x 



27r 



X — z 



-(io"(x) 



Proof. By applying (fT5|) . we can say that if — )■ oo and K — oo, then, 

^ p f poo 

c;TEg(z; LOk) ^ 27r {£-eo)JoiuJk\^ - z\)da{x.) ^ 2n / / (£-£:o)^o(w|x - z|)(ia;(io-(x). 

k=l k=i ^'^ "^0 

Since following infinite integral of the Bessel function formula holds (see formula 11.4.17 

(page 486)]), 

poo 

/ Jnit)dt = 1 (16) 

Jo 

for 9^e(n) > — 1. Then, applying change of variable a;|x — z| = t in (ITB]) for n = yields 

Mt) 



Jo(u;|x — z\)duj 



-dt 



X — z 



X — z 



Thus, dT^s{z]UJk) becomes 



(irE^z; Wfc) 



/ (^-^o)- 



27r 



X — z 



j-(icr(x) 



and by taking the maximum value, we can obtain the desired result. Similarly, dx^fj^^z; Uk) 
can be written as 



2( i-lV,.t(x) + 2f— -4 ld/-n(x) 



j-dcr(x) 



On the basis of Theorems 13.21 and 13.41 we can explore some properties of normalized 
multi-frequency topological derivative imaging function f lTOj) . summarized as follows: 



1. io(x) and cosx reach their maximum value 1 at x = and x = 2n7r, respectively, for 
n = 0, ±1, ±2, ■ ■ ■ . Therefore, E(z; K) plots its maximum value at z, which satisfies 

^idi ■ (x — z) = and ^2d; • (x — z) = 2n7r 

for n = 0,±1,±2,---. This implies that points of magnitude 1 (or close to 1) will 
appear at z = x, i.e., along the unknown supporting curve a. Moreover, since 

sin ax 

lim cos bx — > 0, 

a:-s>oo ax 

E(z; K) plots when z is far away from x. 

2. A(a;; uk) — A(a;; Ui) has properties similar to those of joi^x) cosx, except for less oscil- 
lation. Hence, E(z; K) plots its maximum value at z = x G a. 

3. jo{ax) cos{bx) and A(a;; w/^) — A(x; ui) have their minimum values at two points Xi and 
X2, symmetric with respect to x, refer to Figures [2] and |3l respectively. This implies 
that the map of E(z; K) contains its minimum values in the neighborhood of a so 
that the location of the supporting curve is clearly identified by looking at points of 
maximum and minimum values. 

4. Applying multi-frequency (i.e., K is sufficiently large enough) will guarantee a better 
imaging result than single frequency (i.e., K = Vj. Moreover, it is expected that 
applying a postprocessing operator introduced in |5| will yields a better result. 

5. The map of E(z; K) accurately yields the location of z = x G a when we apply a large 
number of K and L. 



4. Numerical results and discussions 

4.I. General configuration of numerical simulations 

Some numerical simulation results are presented herein. For simplicity, we consider the 
dielectric permittivity contrast case only. The homogeneous domain Q is chosen as a unit 
circle centered at the origin in R^, and three aj specify the thin inclusions as 

(Ti = {(s — 0.2, — 0.5s^ + 0.5) : s G [—0.5, 0.5]} (curve with constant curvature) 

o"2 = {(s + 0.2, + — 0.6) : s G [—0.5, 0.5]} (curve with nonconstant curvature) 

(T3 = {(s,0.5s^ + 0.1sin(37r(s + 0.7))) : s G [-0.7,0.7]} . (oscillating curve) 



-5 -4 -3 -2 -1 1 2 3 4 5 

X-axis 

Figure 2: Graph oi y = jo{ax) cos{bx) for a — 2 and b — 10. 
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Figure 3: Graph of y = A(x; ujk) — A(a;; wi) for ujk = ^ and uji = 



The thickness of the thin inclusion Tj is set to 0.02, and parameters Eq, fio are chosen as 
1. Let Sj and fij for j = 1, 2, 3 denote the permittivity and permeabihty of Tj, respectively. 
The applied frequency is selected a.s Uk = at wavelength X^, k = 1,2, ■ ■ ■ , K and L = 4 
different incident directions 

:= ( cos ^ ^ ^ ,sm ^ ^ ' j, l = l,2,---,L, 

have chosen. In order to show the robustness of the proposed algorithm, a white Gaus- 
sian noise with 15dB signal-to-noise ratio (SNR) added to the unperturbed boundary data 
u^^\'K.;Uk) via a standard MATLAB command 'awgn'. Throughout this section, only both 
permittivity and permeability contrast case is considered, and we select Sj = fij = 5 for 
j = 1, 2, and 3. 
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4-2. Numerical results and discussions 

First, let us consider the influence of the number of frequencies K. For this purpose, we 
choose a thin inclusion Fi and compare maps of E(z; K) for K = 1,5, 10, and 16. From the 
results in Figure HI it is difficult to recognize the shape of Fi when we apply K = 1 or K = 5 
because so many unexpected points of large magnitude are distributed on i7\Fi. However, 
when we apply sufficiently large K, it is easy to recognize the shape of Fi. Based on the 
obtained image, ii' = 16 is a good choice; hence, we will adopt K = 16 different frequencies 
in this section. It is interesting to observe that when K increases, the points of minimum 
value of E(z; K) appear in the neighborhood of Fi. 

Maps of E(z; K) are shown in Figure [5] when the thin inclusion is F2. Similar to the 
imaging of Fi, we can identify F2 when the value K is sufficiently large. 

Let us apply the imaging function to F3 under the same configuration as the above 
examples. Although only four points of F3 are clearly identified, E(z; K) offers an acceptable 
result for an oscillating inclusion by comparing the result in |23|, Figure 5]. 

11 





12 



-1 



-0.5 0.5 1 

X-axis 



-1 



-0.5 0.5 1 

X-axis 



Fi gure 7: Maps of IE(z; for L — 4 with K — 4 (left) and K — 16 (right) when the thin inclusion is Fm- 



One advantage of topological derivative is its straightforward application to the imaging 
of multiple inclusions. Figure [7]shows the map of E(z; K) for imaging multiple thin inclusions 
Tm = Tmi U = Fi U r2 with £i = £2 = 5 and /Ui = /i2 = 5. Unlike to the previous single 
inclusion cases, although the existence of two inclusions can be recognized, it is difficule to 
identify their true shape. 

Figure [S] shows the map of E(z; K) under the same configuration as the previous example, 
except for different material properties, £1 = /ii = 5 and S2 = ^2 = 10. Note that in the 
existence of M— different thin inclusions, Theorem 13.21 becomes 



Ei(z;ir) = ^5^ / (£„-£o)jo( 

L M - 
1=1 m=l 



■ (x — z) ) COS ( ■ (x — z) acr(x) 



2( )d,.t x) + 2 --^ d;-n x) 



X Jo 



1 1 

/^m /^O 



1 /^m 



y- 



d/ ■ (x — z) cos 



/^o ^5 

UJR + 



d/ ■ (x — z) J (icr(x 



Theorems 13 . 31 and 13 .41 can be written in a similar manner. Hence, it is true that if an inclusion 
(here, Fi) has a much smaller value of permittivity or permeability than another (here, F2), 
this inclusion does not significantly affect the scattered field, and as a consequence, the value 
of E(z; K) for z G Fi will be smaller than E(z; K) for z G F2. 

An improvement can be realized by simply making L as large as possible. Figures [9] and 
[To] are maps of E(z; K) for L = 16 in the existence of a single thin inclusion. By comparing 
Figures m and El the shape of F^ appears more accurate than the L = 4 case. Note that if 
one can apply a large number of incident directions L, the number of applied frequencies K 
can be reduced, refer to Figure fTOl Figure [12] shows the map of E(z; K) with K = L = 16 
in the existence of multiple thin inclusions. As expected, good imaging results are obtained. 
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Figure 10: Maps of E(z; K) for L — \Q with K = A (left) and K — L — \Q (right) when the thin inclusion 
is Ta. 



Now, let us compare E(z; K) with two well-known non-iterative algorithms, Multiple 
Signal Classification (MUSIC) and Kirchhoff migrations (see Appendix |A] for corresponding 
algorithms). Figure [TT] shows the imaging result of MUSIC and Kirchhoff migrations for 
L = 16 when the thin inclusion is without noisy data. From Figures dO] and [TTl we can 
observe that because of the small value of Z0, a good result cannot be obtained via MUSIC 
and Kirchhoff migrations, but E(z; K) yields a good result. 

4-3. Producing a good initial guess for applying iterative algorithms 

From the results presented in the previous section, we can generate a good initial guess 
for iterative based reconstruction algorithm ji], [l^, 15, 28|. Borrowing the basic idea of 20 



we assume that the supporting curve aj can be represented as follows 

aj = {zj{s) : s G [aj,bj]} , 
where Zj : [aj, bj] — > is of the form 

^ji^) = I s^^CpTp{s) J , s e [aj,bj]. 
\ p=0 J 

Here, Tp{s) denotes the Chebyshev polynomials of the first kind, defined by the recurrence 
relation 

To{s) = 1 
Tiis) = s 
Tp+i{s) = 2sTp{s) - Tp_i(s). 



-'^If the value of L is sufficiently large enough, good result can be obtained, refer to 21, 25, 13 

15 



0.5 



ns 



I 



I 



1.6 
1.4 
1.2 
1 

0.8 
0.6 
0.4 



0.5 



CO 



-0.5 



I 



I 



700 
600 
500 
400 
300 
200 
100 



-0.5 





X-axis 



0.5 



Figure 11: Imaging result via MUSIC (left) at single frequency uj 
multi- frequency when the thin inclusion is T^. 



and Kirchhoff migration (right) at 





Figure 12: Maps of E(z; K) for K = L = 
permeabilities when the thin inclusion is Fm- 



16 with same (left) and different (right) permittivities and 



16 



From the numerical experience in 20|, Section 7], we use g = 5 polynomials Tp{s), p 



1, 2, ■ ■ ■ ,q, in order to represent aj. The computed coefficients Cp listed in Table [H and the 
corresponding curves cr^"'' are shown in Figure [131 
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0.5235 


-0.2174 
0.1682 


0.0000 
0.2555 


0.0000 
0.0000 


0.0000 
0.0000 


0.3251 


0.4823 


2.6507 



Table 1: Computed coefficients ap of Chebyshev polynomials of the first kind Tp(s), p = 0, 1, ■ 
values of discrete norms Ni(aj), N2(w), and Noo(w) for uj = 



, 5, and 



Let wf!L(x; a;) and ti£)^p(x; u) be the solution of (j2]) in the existence of a true inclusion 
and initial guess F^"'' with supporting curve a^'"*, respectively. Then, due to the difference 
in shape of Tj and F^"", we can define some discrete norms and evaluate them in order to 
investigate the fact that the obtained thin inclusion F^"" is close to the true one Fji for 
x„ e n = 1,2, - ■ ■ , A^, 



L ^ L N 

1=1 1=1 n=l 



^l(xn;c^)-^x(lp(x„;a;) 



L ^ L / N \ i 



2 



/ J 1 1 true V ' / comp V ' / 1 1 ^- y-Ji- } 

1=1 1=1 \n=l 



max |mJ1(x„; u) - ^^^^(x^; u) 



1 ^ 1 ^ 

Notice that in this paper, since is a unit circle, N = 128 different points x^, on the 
boundary dQ are chosen as 

2nn 2nn\ 
cos^,sm^l for n = l,2,--- ,N. 

In Table [H values of Ni(a;), N2(a;), and Noo(<^) = ^ire listed for thin inclusions Fi, 

F2, F3, and Fm- Obtained supporting curves af'^ and corresponding values of discrete norms 
indicate that a good initial guess is obtained and it will be useful for performing complete 
shape reconstruction via an iterative algorithm, for example, level set method introduced in 
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5. Conclusion 

We investigated the applicability of a multi-frequency based topological derivative algo- 
rithm for the imaging of two-dimensional thin, penetrable inclusions embedded in a homo- 

17 





18 



geneous domain. Various numerical results indicate that the proposed algorithm is stable 
even in the existence of random noise, and it is applicable for single and multiple inclusions. 
Although the shapes obtained via imaging results do not yield the complete shape of the 
inclusion with certainty, the iterative reconstruction algorithm can be successfully performed 
by employing them as an initial guess. 

However, the proposed algorithm has some limitations; for example, it cannot be applied 



to the limited-view inverse problems in contrast to the Kirchhoff migration, refer to [6|, |21 



22l . |26|. Therefore, the supplement of a deficiency point of the proposed algorithm will be 
an interesting subject. 

In this contribution, we considered the imaging of thin electromagnetic inclusions when 
measured boundary data is polluted by Gaussian random noise. We believe that the pro- 
posed algorithm can be applied for imaging when the measured data is distorted by random 
scatterers. 



A. MUSIC algorithm and Kirchhoff migration 



In this appendix, we briefly introduce the well-known MUSIC algorithm and Kirchhoff 
mi grat ion. More detailed discussion can be found in various literatures (g], 0, 0, 10, 14, 18, 



21 



22L l25l . l27l | and references therein. 
Same as in section El let Ut]t{'x;u) and ti[2c(x;a;) denote the total and background solu- 
tions of ([2]), respectively. Then, scattered field measured at boundary dfl can be written as 
an asymptotic expansion formula (see [ll| for instance). 



M(2(y; a;) - «lil(y; u) = h [ Vniil(x; u) ■ M(a; x) ■ VAf{^, y; co) 



+ uj'\e - £o)^^lac(x; tj)A/'(x, y 



(17) 



(i(T(x) + o(/i), 



where A/'(x, y; uj) denotes the Neumann function for Helmholtz operator + u'^EqUq in Q 
corresponding to the Dirac delta function — 5(x, y) that satisfies 

VW(x, y; u) + a;2£o/UoAr(x, y; u) = -5(x, y) in n 



du{x.) 



on dVl, 



and a symmetric matrix M(cr; x) is defined as follows: for x G o", let t(x) and n(x) denote 
unit tangent and normal vectors to a at x, respectively. Then 

• M(cr; x) has eigenvectors t(x) and n(x). 

• The eigenvalue corresponding to t(x) is 2 — ^j. 

• The eigenvalue corresponding to n(x) is 2 
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J M 

MO A«o 



Now, let {di}lLi C be a discrete finite set of incident directions and {d^}^^^ C be 
tlie same number of observation directions. Then, applying asymptotic expansion formula 
(fT7|) and performing integration by parts, we can obtain the following normalized boundary 
measurements: 



an \ J (J^\y) 



dS{y) 



h 



an 



VwliUx; u) ■ M{cx; x) ■ VAf{^, y; u 

+ oj'^ie - £o)Mbal(x; w)A/'(x, y; u 
--h 



^ duiil{y;uj) 
^^^^^ du{y) ^'^^^^ 



d(j{x.) 



j {{e - Sq) - di ■ M(a; x) ■ d^^^ e-'"('^^^-'^')-"d(T(x). 



With this, we can generate a Multi-Static Response (MSR) matrix A = (Aj7(x; oj))^^ 



1=1 



C whose element Aji{yi]u) is the collected normalized boundary measurement at obser- 
vation number j for the incident number /: 



"S(y;-)-"S.(y;-))^|^<'S(y), 

an \ J ou[y) 

ioT j,l = 1,2,--- ,L. 

It is worth emphasizing that for a given frequency lo = based on the resolution 
limit, any detail less than one-half of the wavelength cannot be retrieved. Hence, if we 
divide thin inclusion T into M different segments of size of order 4, ffl^V Q^^^. P9j^t ^say , x.^ , 



™ = 1, 2, . . . , at each segment will affect the imaging (.ee Sfla 

For the sake of simplicity, let us set dj 
observation directions configuration, and assume that M < L; then 



10|,|2l|,l22|,l25|,|26|,l27|) 



dj, i.e., we have the same incident and 



A.Kx; co) =hu^ j (^{e - Eq) + di ■ M(a; x) ■ d,^ e'"('^^^+^') "rfa(x) 



\a\ 
M 



M 

E 



(e - eo) + 2 ( - - — ) dj ■ t(x„)d; ■ t(x^) 



+ 2 ( — - ) dj ■ n(xm)di • n(xm) 



^iuj{6.j+di)--x.„ 



where \a\ denotes the length of a. 

Now, let us perform the Singular Value Decomposition (SVD) of A 



M 
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and define a vector w(x; a;) G C"^^^ as 

w(x; u)=(c-il, dl)e^'^'^-^ c ■ (1, d2)e^"'^-^ ■ ■ ■ , c ■ (1, dz.)e*"^-^^ , (19) 

where the selection of c G {0} depends on the shape of the supporting curve cr(x) (see 
271 Section 4.3.1] for a detailed discussion). Then, by defining a projection operator P onto 



the null (or noise) subspace, for L x L identity matrix I^,, 



M 



P(w(x; u)) := II - ^mi^m'^ w)um(x^ 



m=l 



we can construct MUSIC-type imaging functional: 



EMusic(x;a;) = / — (20) 

||P(w(x;a;))|| 



Now, we introduce Kirchhoff migration; 

L 



Ekm(x;w) := |w(x; a;)Aw(x;u;)| = ^ Sm(w)|w(x; a;)u„(x„; a;)||w(x; a;)v„(x„ 



i^J 

m=l 



where w(x;c<;) is defined in (|T9i) . Note that based on the Statistical Hypothesis Testing, 
multi-frequency Kirchhoff migration 



K K L 



EMKM(x;a;fc) := ^EKM(x;a;) = gm(c^fc)|w(x; a;fc)um(x^; a;fc)| |w(x; ujk)vm{^m; ujk] 



k=l k=l m=l 



will yields a more more accurate result than the single frequency case (see 0, H 0, 
for a detailed description). 

B. Proof of Lemma 13.11 

Now, we shall show a proof of Lemma 13.11 For the purpose of simplicity, we set Sq = 
fj,Q = 1, e > Eo and fi > fiQ. 

First, let us explore the structure of dx^si^', i^) in ([ZD- Notice that in this case, e Eq and 
H = fiQ. Since fldj(x; f^) satisfies adjoint problem (|9]), it can be represented by the Neumann 
function A/'(x, y; a;): for z G f2, 

v^liz;u) = j^^ '^^J^M{7.,y-u)dS{y) = j^^ (^n«(y;a;) - ut{T,uj)^M{7.,y,u)dS{y) 
on i^y an ^^^^ 
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Plugging formula (l2Ti) into ([7]) and applying asymptotic expansion formula (1T7|) yields that 
1=1 ^ 



(iTlEe(z; 



1=1 



Wtot(y;^) - Mbac(y; A/'(z,y; w)dS'(y) ]uUc{2;u 



--hu\e - e,)^ty] ( [ [ iziil(x;c^)A/'(x,y;cu)da(x)Ar(z,y;c^)d5(y)klil(z;w^ 
L V Jdn J a J 

--hu?{e - e^)^t ^ ^ (^N(x, z; a;)Miil(x; a;)ni2c(x; o;)^ da(x) 



where 



N(x,z;cj) 



Ar(x,y;w)Ar(z,y;a;)c?5(y). 



(22) 



(23) 



By virtue in |9|, Neumann function A/'(x,y;a;) has a logarithmic singularity. Therefore, it 
can be decomposed into the singular and regular functions; 



A/'(x, y; a;) = -— In |x - y| + 7^(x, y; w), 



(24) 



where 7^(x, y; a;) G C^'" in both x and y for some a with < a < 1 (see jo], [loj for instance). 
Since x G cr and y G there is no blow up of A/'(x, y; w), and it can be bounded by 

|A/'(x,y;a;)| < ;^ln|x-y| + |7^(x,y;a;)| < -J-lndiam(fi) + max |7^(x, y; a;)|, 

where diam(f2) denotes the diameter of f2. Hence, applying Holder's inequality yields 

|N(x,z;a;)| < ( — Indiam(fi) + max|7^(x,y;a;)| ) / \M{^,y\uj)\dS{y). 
V27r / ho. 

From the fact that y G 9f2 and z G ^2, we must consider the singularity of M(T^^y\bS) at 
z = y in order to analyze (1251) . For handling this singularity, for a fixed small constant 
p > 0, generate a ball i?(z, p) of center z and radius p such that 

5(z,p) nr = 0. 

Then by separating the boundary into 5fi = d^s U Sfi^ (see Figure [H]) , where 

afis = finas(z,p) and a^ij = a^]\(^]^as(z,p)). 
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Figure 14: Illustration of dfts (red-colored line) and dil]^ (blue colored line). 



Then 



/ m^,y;io)\dS{y)<^ [ \n\z-y\dS{y)+ [ \n{z,y;co)\dS{y) 
Jan Jan Jon 

<^ lim ( [ In |z - y\dSiy) + [ In |z - y\dS{y) 



+ max \7l{z, y; uj)\\ength.{dQ) 
<— lim ( plnp + (length(c?f2) — p) In |length(9n)| 



■ 2n p^04 
+ max |7^(z, y; a;) |length(9i7) 

=length(a(]) In |length(an)| + max |7^(z, y; a;)| . 

Here, length((9fi) denotes the length of dQ. Therefore, we can say that N(x, z; u) is bounded 
by 

|N(x, z; w)| < Indiam(fi) + max |7^(x, y; a;)|^ x 

length(5f2) ( 77- In |length(9i7) | + max |7^(z, y; | ] < +00, 



271 

and there is no blow up of N(x, z; u). 

In this paper, the background solution is selected as Wbic(x; u) = e'"^*^''^. Hence, (I22|) can 
be written as 

drEeiz; u) = huj'^{e - eo)^t ^ j (n{x, z; u)u^^l{x; uj)u^^l{x; uj)^ da{x) 

! 1 ^ cr 



1=1 
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Applying the same process to the magnetic permeabihty contrast case {e — Eq and 
H 7^ Ho), we can obtain the following structure of dT^i^{z; u): 



L 

1=1 
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